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Abstract 


In this report the problem we are going to study, is the interpolation of a set. 
of points in the plane with the use of control theory. We will discover how 
different systems generate different kinds of splines, cubic and exponential, 
and investigate the effect that the different systems have on the tracking 
problem. Actually we will see that the important parameters will be the two 
eigenvalues of the control matrix. 
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1 Introduction 

I would like to begin by thanking my advisor Pr. Anders Lindquist for 
initiating contact with Texas Tech. I would also like to thank Pr. Clyde F. 
Martin for being my advisor at Texas Tech. 

In this report we will look at a way to store signatures. We want to 
do this by storing only a minimal amount of points on the signature curve, 
and still be able to reconstruct the curve by interpolating these points. The 
interpolation will be performed by splines, and we will look at the common 
splines-problem from the control theory point of view. We can construct a 
trajectory of a system that passes through a specified set of points, and thus 
interpolate the points. 

Two questions that need to be answered arise. First, when is it possible 
for the system to pass through the points? Second, when there are many 
ways to accomplish that, what sort of conditions should we demand that the 
system fulfill in order that we get a unique solution. 

The question of when it is possible to interpolate the points will be an- 
swered in the general case in section 2 Reachability and for our particular 
system in section 3 The System . 

An algorithm to find the solution is developed in section 4 Derivations . 

The choice of boundary conditions is discussed in section 5 Boundary 
conditions. 

In section 6 Results the results of tests done upon parametric curves are 
displayed and discussed. 

A summary in Swedish is provided in section T Resume - Summary in 
Swedish . 

The programs I have been using are included in section 3 Programs. In- 
cluded among the Matlab programs is an altered version of the original Mat- 
lab program quadS. 

When we have answered the two questions, we have to decide what kind 
of system we will use for the interpolation. We can easily imagine that we 
would get to completely different curves if we asked a pedestrian to walk 
through a set of points and if we asked a cyclist to ride his bike through the 
same points. In the first case we would get (if we suppose that man is lazy), 
linear interpolation, and in the second case we would get a smooth rounded 
curve. This is the same as in our case where we have exponential and cubic 
parametric splines. 
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2 Reachability 

In this section we will determine under what circumstances it is possible to 
take a time invariant linear system from a point Xo at time to to a point Xi 
at time 1 1# It is a vital property to us, because, in order to interpolate we 
have to be able to pass through the points. We will call a system completely 
reachable if it has the property that this can be done in any positive time 
for any two points. 

This is a classical control theory question, and it is answered by the 
following well known theorem, which was at least implicitly discovered by 
Kalman. 


Theorem 2.1 Suppose that the system below is given , 

f x = Ax + B u 
l y = Cx 


( 2 . 1 ) 


where A is n X n and B is n x k. Then it is completely reachable iff T = 
[By ABy A 2 By . . . , A n ~ l B], is full rank . 

In order to prove this and understand how the reachability concept can be 
characterized by the matrix T, we will have to look at the general solution 
of equation 2.1. 

x(t) = e^- to >xo + /' e A ^- f) Bu k {s)ds (2.2) 

JtQ 

In order to have the desired state x x at time t x the following equality must 
be satisfied. 

x t = e^“-‘ 0) Xo + T e Ait '- 3) Bu k (s)ds (2.3) 

Jt 0 

The question of reachability is now easily seen to be the question of whether 
there are any solutions to the mapping L \ U ► R n such that 

Lu 4 f tx e A ^~^Bu{s)ds = Xl - e A(t, -‘ o) Xo = d (2.4) 

Jt 0 

Since we recognize L as a linear operator, it is as always very fruitful to use 
a theorem from the general theory of functional analysis [4, p.250]. 
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Theorem 2.2 If X, Y are complete inner-product spaces and A:X Y is 
a linear continuous operator then 

3m A — CJm^i4* 


We know that R n is a Hilbert space, but we have to look at what kind of 
space U is. We choose to introduce the following inner product for U 

(u, v)^ = f 1 u(f)'v(*)di 
Jt o 

and it can be checked that U becomes a Hilbert space. VVe know that there 
is a adjoint operator L m :R n >-*U such that 

(d,Iu)Rn = (Z"d,U)tf 

and we get the adjoint L‘ of the mapping L through this equation 

(d,Iu) R n = d' ^ e Mt '~ ,) Bu(s)d3 

Jto 

= Jt ^B'e A '^ l ~^dj u(s)ds = (L~d,u)u 

We thus have a linear mapping W = LL m : R n t-*R n , that is, it is actually 
only a matrix operator. 

W = LL m = r 

Jto 

With only the basic knowledge about matrices we will now be able to prove 
the following le mm a 

Lemma 2.1 Let A be n x n and B be n x k. Then, for all t 0 , t x such that 
to < t\ we have 

amW(f 0 , U) = Jm [j B , AB, A 2 B , .... A n ~ l B] 

Proof: We will do this by showing that Jmr C 3mW and JmlV C Jmr. 
pmT C 3mW] 
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Let a G £&x W, which implies that 0 = a'W’a = a ^ J ^a ds, 
i.e. 

P \B'e A ' {tl - 3 ) a\ds = 0 

Jto 1 1 


from which it follows that 

BV 4 '( {r,) a = 0 T 


V s € [to, ti] , 


he, 

fi(n-W') J a = 0. 

J= 0 J * 

This implies that [ B , .45, A 2 £, . . . , a = 0. That is, for an arbitrary 

a G £erlV we have a G jSzrF' which implies that JizvW C SizvF' and by a 
theorem from fundamental algebra this equals Cfmr C Jm W . 

pmFF C amr] 

Suppose a G 3vc\XV . Then there exists a x G R n such that a = W x, and 
hence 

a = f^ r 4(<1 - s) J B'e A '^xds 

,to J *o 

from which it is obvious that a G 3m [5, ,45, .4 2 £, . . .] and by Cayley- 
Hamiltons theorem that a G UrnT, which concludes the proof. □ 

The main theorem of this section will follow as a direct consequence of 
the lemma. 

Proof: [Theorem 2.1] By lemma 2.1, Xi — e' 4 ^ 1_f °^Xo — d G R n and 3mr = 
implies complete reachability. □ 
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0 10 0 
0 Ai ai a 2 

0 0 0 1 

0i 0i 0 A o 

This gives us the property 

y = Cx = 

and the system dynamics will look like 

f x = Aii 4- aiy + ct-iij + ui 

= A 2 y + 0ix + 0 2 x + u 2 

And it gives us the following F 

ooio Aj <*2 Ti.r r 13 

1 0 A! a 2 oi202 d" Aj au + Q 2 (Ai 4- A 2 ) r 2>7 F 2 , 8 

oooi 02 a 2 r 3>7 r 3t3 

0 1 02 A 2 01 + 02^1 + A 2 ) oc 202 + A 2 T - 1,7 r 4i8 






where 



4 DERIVATIONS 


8 


Ti.7 = <*202 + ^ i 

Ti.s = QTi + £* 2(^1 + A 2 ) 

r 2,7 — &102 “1" &201 "T O 2 A (2Aj + Ao) + Aj 

r 2,8 = 0^2 + a l(^l + ^ 2 ) + CL 2^: + Of2AiA 2 4- 02^2 

r 3,7 = 01 -f A>(Ai -f a 2 ) 

F 3^3 — 0-202 A o 

I\,,7 = 4* ^1(^1 T ^2) + /?2^1 "T* ,^2^1^2 + ,^2^2 

r 4 t g = Otj/?2 + CX201 -f <*2/^2 (^l + 2Aj) + Aj 

As r is easily seen to have full row rank, by simply looking at the first four 
columns. The class of systems we axe going to consider in this article will all 
have the desired property of complete reachability, by Theorem 2.1. 

4 Derivations 

Given a set of points in the plane {(xcb yo)i ( x i 5 yi)i • • • i ( x niyn)} a^d the 
corresponding time points {t 0 , £ 1 , . . . , t„) we would like to find the control 
functions {u 0 , u t , . . . , u n _i } that takes the system through the points at the 
specified times. 

Let’s study the control : 

As t 6 [<*, 4 + 1 ] the state of the system will be 

x(t) = + f e A{t -’ ] Bu k (s)ds (4.3) 

J tk 

and as we want the state of the system to be x*+i at time £*+ 1 we get the 
following condition. 



x* + i = e' 4( ‘* +,-u) Xfc + f +1 e A{tk +'~ ,) Buk(s)ds 

Jtk. 


(4.9) 
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The solution u* to equation 4.9 that minimizes the norm of the control 
signal is then given by 

Ufe(t) = B'e- A,t (jp* 1 e~ Aa BB'e- A ' 3 ds^ {e~ Atk+l x k+l - 8" i4l ‘x*) (4.10) 

The control would be specified completely by equation 4.10 if we knew the 
whole state-vector at each interpolation point. We know all (xfc, yk) but we 
do not know the ( ik,Vk )♦ To determine the 2 (n + 1) unknowns we have to 
apply some conditions on the solution, and our first choice will be to require 
that the control is continuous, 

Assumption 4.1 

Ufe(£fe+x) = u Jfc+i(tfe+i)» fc = 0...n — 2 


This will give us 2(n - 1) conditions and will leave only four unknowns. We 
will apply the additional conditions on the boundary and we will have to 
come back to this in section 5 Boundary Conditions. 

In manv applications it is just the shape of the trajectory that matters, 
and not the velocity that the system tracks the trajectory with. In these 
cases it makes it much easier to assume that the time between each point is 
a constant. 

Assumption 4.2 Let tk + 1 — tk = h. 

Assumption 4.2 can be used to simplify the integral in equation 4.10. 

/ U+1 e- Aa BB'e~ A ' s ds = {r = s - t k } = 

Jt k 

e~ Atk [ h e~ Ar BB' e~ A ‘ T dr e~ A ’ tk 

^ 

matrixcoTtjtani 

M = f h e- Ar BB'e- A ' T dT 

Jo 


Definition 4.1 
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We can now rewrite expression 4.10 as 

u k (t) = r M x w -Xt) ( 4 . 11 ) 

Using this expression we will now investigate how the continuity condition 
in assumption 4.1 looks. 

u fc (t fe+1 ) = B'e- A ' h M- l (e- AK x k+l - x k ) = 

B'M- l {e- Ah x k +2-Xk+ 1) = u*+i(*i+i) 

We can rewrite this condition using 

Definition 4.2 

Z = M- x t~ AK 
W = e~ A ' h M~ 1 e~ A>l + M~ l 
giving the equation the simple form 


B' {Z’x k - Wx k+l + Zx k+2 ) = 0, * = 0...n-2 (4.12) 

In block diagonal form 


' Z' 

-w 

z ... 

0 0 

0 ' 

0 

Z' 

-W ... 

0 0 

0 

0 

0 

0 ... 

-w z 

0 

_ 0 

0 

0 ... 

Z' -w 

z 


Xo 

Xi 

X 2 


Xn-2 

Xn_l 

X n 


= 0 


(4.13) 


Now, we have to look at what the unknowns axe. The vector in equation 4.13 
is made up of subvectors 

[ Xk " 

_ ik 

yk 

. iik . 



4 DERIVATIONS 


11 


consisting of two knowns and two unknowns. By partitioning the submatrixes 
W, Z as 



: : ; : 


! ! ! * 


’ . . . Z\. . . . ' 

w = 

ltf .1 W.2 W. 3 W.4 

,Z = 

2-1 2.2 2.3 Z.4 

= 

. • . Z'y. ... 
... 23 . ... 


: : : : 


: : : : 


. - • • 2 4 . . . . . 


and using the notations given in the following the definition, we can keep 
the unknowns on the left side and move the position coordinates over to the 
right hand side. 

Definition 4.3 


x ?OJ 


X 


vet 


w? 

wf 


7 b 

"f u 


rB 


7 b 
4 'll 

7 b 

& ri 


X 

y 


X 

. y m 

B' [ 10. 2 

B' [ u».i 

B' [ z . 2 

B' [ 2.1 


W. 4 
U7.3 


a 1 2 2 

W 24 

W 42 

W44 

W 21 

W 2 Z 

W 41 

U>43 


2-4 

Z.3 


-22 ^24 

-42 -44 

-21 z 23 

241 ^43 


B' 


B' 


Zi- 

/ 

Z 22 

242 

2 4 . 


Z 2A 

244 

,,1 

/ 

z 12 

232 

J 


Z 14 

234 _ 
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And we get the system 




0 

0 


-W B 

7 s 
^ u 


0 

0 


za 

-W, B 


0 

0 


0 

0 


0 

0 


-W B Z,Z 0 

-W B Z B 


Zi f 


*0 

v-yc/ 

A 1 

^rVcl 

Ao 


X n-2 

V vei 

X n-1 

v uei 








r v ?os i 
x o 

nT° ' 

1 

-W B 

Z% 

Z r i ... 

-W B ... 

0 

0 

0 

0 

0 

0 


X 1 

X 2 

0 

0 

0 

-W B 

N .. 

0 


X n— 2 

0 

0 

0 

zs 

-\V B 

Zl 


x n~ 1 







•i 

X 


As we evaluate the right hand side, we get a constant vector. Depending 
on what kind of boundary conditions we choose to use, the derivations differ 
from here. We will deal with the most common cases, each case in turn, 
beginning in the next chapter. 


5 Boundary Conditions 

In order to get a unique solution to our problem, we had to apply the con- 
tinuity condition and the boundary conditions. The continuity condition is 
a rather natural condition, but the boundary conditions have to be studied 
more extensively. The four most common choices of boundary conditions in 
the one dimensional case according to [2] are 

1. Zero velocity at the first and at the last point. 

2. Specified starting and ending direction. 

3. Natural boundary conditions, y ft = 0. 

4. Periodical conditions. 
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We will look at how the two dimensional equivalent of choice number 1 
and 2 affect the curves, and for which the same derivation is valid. 

Item 3:s two dimensional equivalent, x = y = 0, requires a derivation and 
will be studied in subsection 5.2. 

We will avoid dealing with condition 4 as it only complicates the calcula- 
tions, and would only be natural and interesting if we were to write the same 
word twice, connected with itself as RogerRoger. 

Because the effect of the boundary conditions are similar at both bound- 
aries we will restrict ourselves to only talking about the starting point. 

5.1 Known velocities at boundary 

We have assumed that we also know and x" e< . Moving these over to the 
right hand side, we get a block diagonal matrix system to solve. 


' -w, B zg ... o o' 

Zg -W B ... 0 0 


Y ue/ 

A 1 

Ao 


fii 

n 2 

0 0 ... -w t B zg 


v ue/ 
X n- 2 


2 

0 0 ... zg -wg _ 


Y vel 

_ x n— 1 . 


f^n-1 


Where 

fJi = -Zjxr + W/xT - Z® XT - W 
fi, = -z3x£, + w - 

o — yS pos , U/B pos rvB poa yB vei 

1 “ -< 6 r /X n _ 2 + VV r X n _ t — ^ rii X n - Zr /u X n 

This can easily be solved, and with a linear increase of time. Having solved 
the system above we now know all the states of the system in the interpolation 
points. We will now use equation 4.11 for the control, and insert it into 
equation 4.8 to get the trajectory. 

x(t) = t A ^U k + f Q ~ t \- AT BB'e- A ' T drM-\e- AK x k ^-y.S) 

for k — 0 ... n — 1 


(5.15) 
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As we can see from equation 5.15 the fundamental matrix and the integral 
is the same for all k and only has to be evaluated for t — t , t between 0 and h, 
once and for all. 

5.1.1 Zero initial velocity 
We can choose to set 
Assumption 5.1 


(~o,yo) — 0 

0=n,yn) = 0 


With this condition, we will let the system start up in whatever direction 
that minimizes the energv norm of the control signal and takes the svstem 
to the second point. 

As we know from one dimensional control theory, a system with a transfer 
function with zeros in the numerator will start off in the opposite direction 
to where it is going. Such undesired properties should certainly be avoided in 
our tracking problem. In the case where a\ = a 2 = = 3 2 = 0, the states 

x and y are independent, yielding two one dimensional transfer functions. It 
is easily seen to have no zeros, which is good. 

Otherwise, we will have to look at the two dimensional transfer function 
given below: 


T(s) 


1 

'S 4 ““ (Ax + A 2 )s 3 T (A1A2 — *r Otll3i)s — Otipi 


X 


s(s — A 2 ) a 1 + ci25 
01 + 02 3 «s($ A| ) 


(5.16) 


This is a bit more tricky, and we will have to find the Q and D of least degree 
that is a solution to the equation 

T(s) = Q{s)D{s)- 1 (5.17) 


and satisfies 


X{s)Q(s) + Y(s)D{s)=I 


(5.18) 
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It i3 easy to verify that the choice 


Q(s) = /, D(s) 


s(s - A x ) -(ai+a 2 s) 
—{/3i + Si*) s{s — A 2 ) 


is a solution to equation 5.17, and there exists X(s) and F(s) so that equa- 
tion 5.18 is satisfied. From Q(s) we can see that the system has no zeros. 
The zero initial conditions should thus not give us any problems. 


5.1.2 Derivative approximation 

Instead of setting the velocity equal to zero, another alternative is of course to 
specify a starting velocity. However, this requires that we make a good choice, 
to avoid situations as exemplified below. Using a bad direction and a high 
velocity boundary condition on y = x 2 , we get the graph of figure 1. Even 
as we are using n=40 to reproduce the curve, the bad boundary conditions 
are still ruining the tracking. 



Figure 1: Trajectory of system tracking y = x 2 , with badly specified starting 
direction 

As discussed in [3, p.86], the fact is that when we set the boundary 
conditions in the parametric case, we do not only specify a direction, but 
also the speed in that direction. The greater the speed, the greater impact 
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the boundary conditions will have on the solution. We are thus forced to 
make a good choice, and we would like the choice to get better the more 
interpolation points we are using. A simple choice that satisfies this is 

Assumption 5.2 

{ioillo) 

(in,Vn) 


which imply that we will set off from the first point in the direction of the 
second point, and arrive at the last point in the direction from the next last 
point. 

Another way of deciding the initial velocity would be to use the same 
technic as in Bezier curves and choose the settings graphically. This would 
probably be the best way to get the desired properties of the signatures. As 
discussed in [3] the choice of boundary conditions will affect the whole curve, 
and the solving of the blockdiagonai system must be done over from scratch, 
making this method a bit slow. If we are going to do this only once to store 
a signature it does not matter. What matters with this method is that it 
adds four more parameters to be stored, and we could probably get equally 
good results just by increasing the number of interpolation points by two. 

5*2 Constant velocity 

Suppose we want to use the boundary conditions 
Assumption 5.3 

(x 0 ,t7o) = 0 

(•^n > j)n ) ^ 0 

This will let the initial direction and constant velocity of the system be 
decided so that the control energy is minimized. Using the system dynamics 
equations 3.7, we get the equation system below. 


Xi - Xo 

h 

X n - X n -1 

h 
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ft : + W, a - U,I) - Zg x£, = (5.23) 

= (u n -w r s - £ “ i jxr+^W-, 

Adding these two equations to equation 4.14 yields a block diagonal system 
to solve. This system is two blocks bigger than the one in section 5.1, but it 
can also be solved with a linear increase of time. And once it is solved, we 
can still use equation 5.15. 

A comparison between the three boundary conditions, BC=1 zero initial 
velocity, BC=2 derivative approximation, and BC=3 zero initial acceleration 
is made in section 6. 
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6 Results 

The first tests with the algorithm were done with matrices A with both 
eigenvalues equal to zero, Aj = 0, A 2 — 0, and c*i = &2 = $2 0- 

This produces cubic parametric splines, and makes the calculation of the 
fundamental matrix easy. The cubic splines produce smooth curves, but were 
also able to reconstruct a cusp much better than you would have guessed at 
first, as shown in figure 2. 



Figure 2: Reconstruction of y = x§ with cubic parametric splines where 
n=10. 


If we look at the function y = we know that this function describes a 
cusp. But if we parameterized it like x = t 3 , y = t 2 we see that both x and 
y axe smooth cubic functions of t , so it is not very remarkable that it can be 
reconstructed well. 

When we used A-matrices with nonzero eigenvalues, and decoupled x 
and y coordinates, i.e. c*i = a 2 = /?i = ft = 0, we were able to generate 
exponential parametric splines with the basis functions l,Aif,e ,6 1 and 

l,A 2 <,e' X2t ,e -A2£ for the x and y coordinates. 

The result of taking big eigenvalues is almost linear interpolation, which 
can be good for certain applications, but not if it is to be used for storing 
signatures. It is quite obvious that the roundness of a persons signature is 
one important factor of it’s characteristics. Therefore, it’s vital that one of 
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data stored on the signature is the eigenvalues of the A-matrix, as can be 
shown in figure 3. 

Figure 4 shows the corresponding control signals. In the Ai = 1, = 1 

case the linear part of the control is dominating, and in the Ax = 100, A 2 = 
100 case the exponential part is dominating. It is also evident that the 
magnitude of the control signals is greater in the case of larger eigenvalues, 
but that is not a problem for our fictitious system. 

In the case of nonzero ax, a?, 3n we get a coupling between the x 
and y coordinates. This will give us very complicated basis functions, like 
polynomials times exponentials times sine- and cosine-functions. 

As with all approximation methods one should always investigate how big 
the errors are. To do this we had to somehow determine the distance between 
the original curve and the interpolating one. The tests were performed on 
known parametric curves, so we had an explicit expression for the points on 
the original curve. We had to try to find the nearest point on the original 
curve to the point on the trajectory. This was solved numerically with the 
w Golden Intersection algorithm' 1 . For the method to work we nave to assume 
two things : 

1. The section of the original curve between the two interpolation points 
nearest in time is convex. 

2. The point on the original curve nearest the point on the interpolating 
curve lies on the section in item 1. 

As an error estimate I have calculated the distance between a number of 
points on the reproduced curve and the original curve and divided with the 
number of points. We have applied this error estimate method on four dif- 
ferent curves and with different number of interpolation points and 40 points 
between each of these. 


points 

BC = 1 

BC = 2 

BC = 3 

n = 10 

6671.1 

5821.1 

4515.9 

n = 20 

102S.4 

765.6 

502.1 

n = 100 

8.8 

5.9 

3.5 


Table 1: fi units of error for unit circle. 
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Figure 5: Graph, of circle, reproduced with n=10 and A x — 10, Ao — 10, 
BC=1,2,3 


6.1 Four test curves 

The first curve we tested was the unit circle. This very round curve was 
tracked best by the cubic splines, but the exponential splines did a good job 
too, as can be seen in figure 5 and table 1. We can also see that the error 
was smallest in the case of zero initial acceleration boundary conditions. 




Figure 6: Graph of y = |x|, reproduced with n=10, n=20 and A x — 10, A 2 — 
10?BC=2 

Next, we looked at a curve with the opposite properties, linear and non- 
differentiable, y = |x|. The error is mainly located between the two points 
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points 

BC = 1 

BC = 2 

BC = 3 

n = 10 

3450.3 

3446.8 

3450.3 

n = 20 i 

972.6 

972.6 

972.6 

3 

II 

O 

CD 

40.7 

40.7 

40.7 


Table 2: n units of error for y = |x|. 


next to the non-differentiable point. As could be guessed the tracking of 
the curve y = \x\ got better the bigger the eigenvalues of the A matrix 
were, the error was reduced to 142 /x units for A* = Aj = 500 in the case of 
BC = 1 and n = 10. but for bigger values the numerical calculations failed. 
Using negative eigenvalues gives the same error as the positive, which can be 
expected since we have a symmetric curve and time interval and by looking 
at the basis functions. The results with different boundary conditions were 
very much alike, as seen in table 2. 

This was the only case were BC 3 did not give us the smallest error. 

Can we get a smaller error with any choice of the coupling parameters. 
Yes, for example by choosing ct\ = — o-i — Pi = ~Pi = we o e ^ the 
error 3249 /x units for n=10. Choosing these parameters could thus be a way 
to reduce the error, but by using n = 12 instead, we got the error 2506 /x 
units. So we do not get a more efficient way to store it, unless we can find 
parameters so we get below 2506 /x units. 


points 

BC = 1 

BC = 2 

BC = 3 

n = 10 

6506.6 

5313.7 

3842.9 

n = 20 

1014.4 

725.4 

460.9 

n = 100 

8.3 

5.3 

3.4 

Table 3: 

pt units of error for 

cycloid. 

points i 

BC = 1 

BC = 2 

BC = 3 

n = 10 

13897.6 

11664.3 

8867.6 

n = 20 

2100.5 

1531.0 

1001.7 

n = 100 

17.7 

11.8 

7.0 


Table 4: /x units of error for prolate cycloid. 
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Finally we look at a cycloid. 

J x = wt — sin xt 
[ y = 1 — cos :ri 

and a prolate cycloid. 

j i = ~t — 2 sin :rt 
[ y = 1 - 2 cos ~t 

It is evident that the cusp and the crossover do not cause any problem, as 
could be expected since we are using a parameterized interpolant. 

We can compare the difference of the graphs in figure 7 and figure S, 
where we are using one crude approximation with n=10 and one extensive 
approximation with n=100. This time it is evident that the error is m ainl y 
located at the boundaries. In table 3 and table 4 we can see that the best 
results came from using constant velocity boundary conditions. 

Looking at table 3 and table 4 again, we see that the error decreases at 
an approximately cubic rate as the number of interpolation points increase, 
which is much better than the quadratic decrease that can be seen in table 2. 

My guess is that this behavior comes from the fact that the curve y = jrj 
does not have a differentiable parameterization. 

6.2 Applied on a signature 

Included as figure 9 is a scanned picture of my own signature. I have tried 
to pick some roughly equidistant points on the signature, (According to the 
scale indicated on the axis.) and used the interpolation algorithm we have 
been studying to reproduce it. The reproduction is made with n = 74, 
Ax = At = 10 and no coupling between the two coordinates. For boundary 
conditions I have chosen to use constant velocity, since it has been the most 
successful condition. 

As can be seen we get a very close resemblance between the original and 
the reproduction. How close is hard to say because we do not have the 
signature given as a parameterized function, therefore we are not able to 
calculate the error as before. 

The things characterizing the signatures, are also the things that are hard 
to recover with the interpolation. Such as the turnover in the connection from 
the "P”, and the cusps in the “r”. To get a good reproduction, an equidistant 
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distribution of the Interpolation points is not enough, more points has to be 
concentrated around the characterizing areas. 


7 Resume - Summary in Swedish 

Lagring av signaturer kan goras pa manga satt. \ i har valt att lagra ett antal 
punkter pa signaturen, och sedan reproducera denna genom att interpolera 
punkterna med splines. 

Genom att anvanda valkanda resuitat inom systemteori sa kan man gener- 
era olika sorters splines genom att andra pa nagra parametrar. Jag har nvtt- 
jat denna metod for att generera parametriserade splines i planet. Man inser 
snart att man maste infora randvillkor pa losningen. och valet av dessa far 
inte ske hur som heist eftersom de paverkar hela losningen. Darfbr har jag 
lagt ner en hel del arbete just pa denna punkt. De basta resulteten har jag 
erhallit genom valet att ha konstant hastighet vid andpunktema. 

En av de saker som karaktariserar handstilar ar dess rundhet. Denna kan 
ges en direkt oversattning i egenvardena till systemmatrisen, och vi kommer 
alitsa valja att lagra dessa utover punkterna pa signaturen. 
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8 Programs 

The programs in Matlab and Maple that were used to implement the algo- 
rithm developed in this report follow. 

8.1 Matlab Program 

To make it easier to understand the structure of the program, the following 
flow charts describe how the programs are traversed. 



The loops marked with an unfilled circle is only available when the inter- 
polated points are given by the parameterized function (XX(t),YY(t)). 
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function error=run(BC,nn) ; 

’/.BC ’/, Type of boundary conditions 

’/, 1 = zero initial and final velocity. 

'/, 2 = beading for first point, last. 

’/, 3 = zero initial and final acceleration. 

'/, If nn is specified, runs alg with nn points given by XX, YY. 
'/, Otherwise runs alg with points given by ginput. 
global n h alphal alpha2 betal beta2 lambdal lambda2 
global errorcalc ctrlsignal 
dg 


Setting of parameters 
alphal =0; 
alpha2=0 ; 
betal=0; 
beta2=0 ; 
lambdal=10 ; 
lambda2=10 ; 

*/, Decides what steps are going to be made 
errorcalc=l; '/, error estimation 

ctrlsignal=0; V. plotting of control signals 

if nargin == 1 
[x,y] =ginput ; 

R ( 1 , : ) =x ’ ; R(2 , : )=y ' ; 
else 

R=points(nn) ; 
end; 

n=length(R)-l ; */. Number of interpolationpoints -1. 

h=2/n; '/, Time inbetween points 

Plots a circle at all the points thats interpolated 
hold on 
for i=l:n+l 
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plot (R(l , i) ,R(2,i),’o') 
end; 

7.7, Calling alg with the prepared data VI. 
if BC == 3 

error=alg2(R) ; 
else 

if BC == 2, xvelO=(R( : ,2)-R( : , 1) )/h; 

else xvel0=zeros(2,l) ; 
end; 

if BC == 2, xveln=(R( : ,n+l)-R( : ,n))/h; 

else xveln=zeros (2 , 1) ; 
end; 

error=alg(R,xvelO, xveln) ; 
end; 

7.7. Loop to allow graphic alteration of BC. 7.7. 

b=input(’"l" for graphic mod of BC, "0" to quit 
while b == 1 , 

xvel0=10*(ginput(l) ’-R( : , 1)) ; 
xveln=-10*(ginput(l) ,- R( : ,n+l)) ; 
dg 

hold on 
for i=l:n+l 

plotCR(l.i) ,R(2,i) , ’o’) 
end; 

error=alg(R,xvelO,xveln) ; 

b=input( , "l" for graphic mod of BC, "0" to quit ’); 
end; 

end; 


function R=points (nn) ; 

% Forms R with the help of XX, YY. 
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7. R = 2* (nn+1) -matrix, 
global a b h 


a=l ; b=a; 
h=2/nn; 

for i= 0:nn 

R( : ,i+l)=[XX(-l+i*h) ; YY(-l+i*h)] ; 
end 
end; 


function res=XX(t) ; 
global a b 

res= a*t*pi-b*sin(t*pi) ; 
end; 


function res=YY(t) 
global a b 

res= a-b*cos (t*pi) ; 
end; 


function error=alg(R,xvelO .xveln) ; 

7. R = Matrix of interpolationpoints , xO ... xn. 

7. first row = x-coordinates , second row = y-coordinates . 

7. xvelO, xveln = Boundary conditions 

global a b A B C n h alphal alpha2 betal beta2 lambdal lambda2 
global errorcalc ctrlsignal 

7.7. The System 7.7. 

A=[[ 0 1 0 0] 

[ 0 lambdal alphal alphal] 
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[ 0 0 0 1 ] 

[ betal beta2 0 lambda2]] ; 

B= [ [0 0] 

[1 0 ] 

CO 0] 

CO 1]]; 

C=[[l 0 0 0] 

[0010]]; 

7,7, Calculation of the integral from 0 to h in m steps 7.7. 
m =40; 7, number of points between interpolationpoints 

tol=le-08; 7. the numeric error tolerance 

Mtau( : ,l:4)=zeros(4) ; 
tau=0 ; 
for j=l:m 
oldtau=tau; 
tau=oldtau+h/ m ; 

Mtau(: ,4*j+l : 4*j+4)= quad8mod( ’ int ’ , oldtau.tau.tol) 

+ MtauC : ,4*j~3:4*j) ; 

end; 

M=Mtau( : ,4*m+l :4*m+4) ; 


7.7. Forming of the Matrixes for the Blockdiagonal system 7.7. 


e_Ah=expm(-A*h) ; 
Minv=inv(M) ; 
ZZ=Minv*e_Ah; 
WW=e_Ah ’ *ZZ+Minv ; 


WL=[WW(2,2) WWC2.4); WW(4,2) WW(4,4)] ; 7. Partitioning matrixes 

ZLU=[ZZ(2,2) ZZ(2,4); ZZ(4,2) ZZ(4,4)] ; 

ZLL= [ZZ(2 , 2) ZZC4.2); ZZ(2,4) ZZC4.4)] ; 

WR= [WW (2,1) WWC2.3); WW(4,1) WW(4,3)] ; 

ZRU= [ZZ(2, 1) ZZC2.3); ZZ(4,1) ZZ(4,3)] ; 

ZRL=[ZZ(1 ,2) ZZ(3,2); ZZ(1,4) ZZC3.4)] ; 
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7.7. The boundary conditions 7.7. 

xvel( : , l)=xvelO; 
xvel( : ,n+l)=xveln; 

7.7. Forming of the right side of the Blockdiagonal system 7.7. 
for i=2:n 

Omega(:,i)=ZRL*R(: ,i-l)-WR*R(: ,i)+ZRU*R(: ,i+l) ; 
end; 

Omega ( : ,2)=0mega( : ,2)+ZLL*xvel( : , 1) ; 

Omega( : ,n)=0inega( : ,n) +ZLU*xvel( : ,n+l) ; 

7,7 Gausselimination to produce upper triangular system 7.7. 

DD ( : ,3:4) =WL ; 
for i=3:n 

zd=ZLL*inv(DD( : ,2*i-3 : 2*i-2) ) ; 

DDC: ,2*i-l:2*i)«WL-zd*ZLU; 

Qmega(: ,i)=0mega( : , i)+zd*0mega( : ,i-l) ; 
end 


7.7. Backsubstitution to solve for the xvel 

xvel(: ,n)=DD(: ,2*n-l : 2*n) \0mega( : ,n) ; 
for i=n-l:-l:2 

xvel( : ,i)=DD( : ,2*i-l:2*i)\(ZLU*xvel(: , i+l)+0mega( : ,i)); 
end; 

7.7. Making of state vectors 7.7. 
for i=0:n 

x(: ,i+l)=[[R(l,i+l)] 

[xvel(l,i+l)] 

[R(2,i+1)] 

[xvel(2,i+l)]] ; 
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function error=alg2(R) ; 

7, R = Matrix of interpolationpoints , xO ... xn. 

’/, first row = x-coordinates , second row = y-coordinates . 

7. BC = acceleration in x and y direction are both = 0. 

global a b A B C n h alpbal alpha2 betal beta2 lambdal lambda! 
global errorcalc ctrlsignal 

7.7. The System 7.7. 

A=[[ 0 1 0 0] 

[ 0 lambdal alphal alpha2] 

[ 0 0 0 1 ] 

[ betal beta2 0 lambda2]]; 

B= [ [0 0] 

Cl o] 

[0 0] 

CO 1]]; 

C=[[l 0 0 0] 

[0010]]; 

Calculation of the integral from 0 to h in m steps 7.7. 
m=40; 7. number of points between interpolationpoints 

tol=le-08; 7. the numeric error tolerance 

Mtau( : , 1 :4)=zeros(4) ; 
tau=0; 
for j=l:m 
oldtau=tau; 
tau=oldtau+h/m; 

Mtau ( : , 4*j + l : 4* j +4) = quad8mod( ' int ’ , oldtau , tau , tol) 

+ Mtau( : ,4*j-3:4*j) ; 

end; 

M=Mtau( : ,4*m+l :4*m+4) ; 
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7.7. Forming of the Matrixes for the Blockdiagonal system 7.7. 

e_Ah=expm(-A*h) ; 

Minv=inv(M) ; 

ZZ=Minv*e_A_h; 

WW=e_Ah ’ *ZZ+Minv ; 


WL= [WWC2.2) WW(2,4) ; WW(4,2) WV(4,4)]; 7. Partitioning matrixes 

ZLU= [ZZ(2 , 2) ZZ(2 , 4) ; ZZ(4,2) ZZ(4,4)] ; 

ZLL=[ZZ(2,2) ZZ(4,2) ; ZZ(2,4) ZZ(4,4)]; 

WR= [WW (2,1) WWC2.3); WW(4,1) WW(4,3)] ; 

ZP.U= [ZZ(2, 1) ZZ(2,3) ; ZZ(4,1) ZZ(4,3)] ; 

ZRL-CZZC1.2) ZZ (3,2) ; ZZ(1,4) ZZC3.4)] ; 

Ulu= [Minv(2 , 2) Minv(2,4) ;Minv(4,2) Minv(4,4)]; 

Uru= [Minv(2 , 1) Minv(4, 1) ;Minv(2,3) Minv(4,3)]; 

Vl= [lambdal alpha2 ; beta2 lambda2] ; 

Vr=[0 alpha 1 ; betal 0]; 

7.7. Forming of the right side of the Blockdiagonal system 7.7. 
for i=2:n 

0mega(: ,i)=ZRL*R(: ,i-l)-WR*R( : ,i)+ZRU*R(: ,i+l); 
end; 

OmegaC : , l)=(Vr-Uru)*R( : , 1) + ZRU*R( : ,2) ; 

OmegaC : ,n+l)=ZRL*R( : ,n) - (WR-Uru+Vr)*R( : ,n+l) ; 

7 . 7 . Gausselimination to produce upper triangular system 7 . 7 . 

DD(: ,l:2)*Ulu-Vl; 
for i=2:n+l 

zd=ZLL*inv(DD( ; ,2*i-3:2*i-2)) ; 

DD(: ,2*i-l:2*i)=WL-zd*ZLU; 

OmegaC : , i)=0mega( ; , i)+zd*0mega( : , i-1) ; 
end 

DD ( : , 2*n+l ; 2*n+2) = (WL-Ulu+Vl) - zd^ZLU; 
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Backsubstitution to solve for the xvel 

xvel( : ,n+l)=DD( : ,2*n+l : 2*n+2) \0mega( : ,n+l) ; 
for i=n:-l:l 

xvel( : , i)=DD( : ,2*i-l : 2*i) \ (ZLU*xvel( : , i+l)+C:nega( : ,i) ) ; 
end; 

7.7. Making of state vectors 7.7. 
for i=0:n 

x(: ,i+l)=[[R(l,i+l)] 

[xvelCl ,i+l)] 

CR(2,i+l)] 

[xvel (2, i+1)] ] ; 

end; 


7.7. Plotting of trajectory 

7.7. and error estimate calculation 

sumnonn=0 ; 
hold on 
for j=0:m 

eAtau=expm(A*j*h/m) ; 
for i=0:n-l 

entry=eAtau*(x( : ,i+l)+Mtau( : ,4*j+l:4*j+4)*Minv* 

(e_Ah*x( : , i+2)-x( : , i+1) ) ) ; 
plot (entry (1) , entry (3 ) , ’ 
if errorcalc 

sumnorm = sumnorm + dist(entryCl) ,entry(3) ,i,h) ; 
end; 
end; 
end; 

7.7. Plotting of the control signals 7.7. 


if ctrlsignal 
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pause, clg 

subplot (2 , 1 , 1) .hold on 
subplot (2 , 1 , 2) .bold on 
for i=0:n-l 
for j=0:m 

csignvec ( : , j+1) =3 ’ *expm(-A ' *j *h/ 40) *Minv* 
(e_Ah*x( : , i+2)-x( : , i+1) ) ; 

end; 

subplot (2, 1,1) ,plot(i*h:h/m: (i+1) *h, csignvec Cl , :)) 
subplot (2, 1,2) ,plot(i*h:h/m: (i+l)*h, csignvec (2, :)) 
end; 
end; 

err o r= s umno rm/n/m ; 
end; 


function [Q,cnt] = quad8mod(F,a,b,tol) 

‘/.Alteration of the original matlab toolbox program. QUAD8 
7, Numerical evaluation of an integral, higher order method. 

7, Q = QUAD8 ( ’ F ’ , A , B , TOL) approximates the integral of F(X) 

7. from A to B to within a relative error of TOL. 

7. ’F J is a string containing the name of the function. 

7, The function must return a 4*4-matrix output value if 
7. given an input value . 

7. Q = Inf is returned if an excessive recursion level is 
7. reached, indicating a possibly singular integral. 

7. QUAD8 uses an adaptive recursive Newton Cotes 8 panel rule. 

7. Cleve Moler , 5-08-88. 

7. Copyright (c) 1984-94 by The MathWorks , Inc. 

7. [Q,cnt] = quad8(F,a,b,tol) also returns a function 
7, evaluation count. 

7. Top level initialization, Newton-Cotes weights 

w= [3956 23552 -3712 41984 -18160 41984 -3712 23552 3956]/14175; 

x=a + (0:8)*(b-a)/8; 
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7, set up function call 
for i=x 

y = [y fevalCF, i)] ; 
end; 

7. Adaptive, recursive Newton-Cotes 8 panel quadrature 
QO = zeros (4) ; 

[Q.cnt] = quad8stpmod(F, a.b , tol , 0 , w.x.y , QO) ; 
cnt = cnt + 9; 
end ; 


function [Q.cnt] = quadSstpmodCF.a.b.tol.iev.w.xO.fO.QO) 
7,Alteration of the original mat lab toolbox program. QUAD8ST? 

7. Recursive function used by QUAD8 . 

7. [Q.cnt] = quad8stp(F,a,b,tol,lev,w,f ,Q0) tries to approximate 
7. the integral of f(x) from a to b to within a relative error 
7. of tol. 

7. F is a string containing the name of f. The remaining 
'/, ar guments are generated by quad8mod or by the recursion. 

7. lev is the recursion level. 

7. w is the weights in the 8 panel Newton Cotes formula. 

7. xO is a vector of 9 equally spaced abscissa is the interval. 

7. fO is a matrix of the 9 function values at x. 

7. QO is an approximate value of the integral. 

Cleve Moler, 5-08-88. 

7, Copyright (c) 1984-94 by The MathWorks, Inc. 

LZ7MAX = 10; 

7. Evaluate function at midpoints 
7. of left and right half intervals, 
x = zeros (1 , 17) ; 
xCl : 2: 17) = xO; 

x(2:2: 16) = (x0(l:8) + x0(2:9))/2; 


f ( : ,1:4)= f0( : , 1 : 4) ; 
for i=l : 8 
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f ( : ,8*i-3:8*i) = f eval(F,x(2*i)) ; 
f ( : ,8*i+l :8*i+4) = f 0 ( : , 4*i+l : 4*i+4) ; 
end; 

7. Integrate over half intervals, 
h = (b-a)/16; 

Q1=0 ; Q2=0 ; 
for i=l:9 

Q 1 = qi + h*v(i)*f ( : ,4*i-3:4*i) ; 

Q2 = Q2 + h*v(10-i)*f ( : ,69-i*4:72-i*4) ; 
end; 

Q = Ql + Q2; 

7. Recursively refine approximations, 
if norm(Q - QO) > tol*norm(Q) & lev <= LEVMAX 
c = (a+b)/2; 

[Ql.cntl] = 

quad8stpmod(F, a, c , to 1/2 , lev+1 ,w,x(l:9),f(: ,1:36) ,Q1) ; 
[Q2,cnt2] = 

quad8stpmod(F, c ,b , to 1/2 , lev+1 , v ,x(9 : 17) ,f ( : ,33:68) , Q2) , 
Q = Ql + Q2; 

cnt = cnt + cntl + cnt2; 

end 
end ; 


function res = integrand(v) 
global ABC 

e_AvB=expm(-A*v) *B ; 
res = e_AvB*e_AvB’ ; 
end; 


function [d] =dist(xx,yy , i ,h) ; 

7. Initiating search algorithm. 
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8.2 Maple Program 

with(linalg) ; 
with (student) ; 

alpha! :=1; 
alpha2 : =1 ; 
betal : =0 ; 
beta2 : =0 ; 
lambdal : =100 ; 
lambda2 : =100 ; 

a: =1 ; 
b : =a; 
h: =0 . 2 ; 
n: =22 ; 
m: =15 ; 

R: =vector(n+l) ; . 

XX:=t->a*t*Pi-b*sin(t*Pi) ; 

YY : =t->a-b*cos (t*Pi) ; 
for i from 0 to n 
do 

R[i+1] :=matrix( [[XX(-l+i*h)] , 

[YY(-l+i*h)] ] ) ; 
od; 

A:=matrix([[ 0, 1, 0, 0], 

[ 0 , lambdal , alphal , alpha2] , 
[ 0 , 0 , 0 , 1] , 

[ betal, beta2, 0, lambda2]]); 
B : =matrix( [ [0 , 0] , [1 ,0] , [0 ,0] , [0 , 1] ] ) ; 

C : =matrix( C [1 , 0,0 ,0] , [0 ,0 , 1 ,0] ] ) ; 

aliasCId = &*0) 

Aprim : =transpose (A) ; 

Bprim: =transpose(B) ; 
e_At :=t->exponential(-A*t) ; 
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e_Aprimt :=t->exponential(-Aprim*t) ; 
eAt : =t->exponential(A*t) ; 
e_.Ah: =e_At(h) ; 


integranden: = proc (v) 

evalm(e_At(v) k* B k* Bprim k* e_Apriot(v)) ; 

end ; 

integrera:= proc (funk) 
global 7 ; 
int(funk,v) ; 
end; 

map ( integrera , integranden (v) ) ; 
int egralen : =map ( s implif 7 , " ) ; 

evaluera:=proc (funk) 
global tau; 
subs (v=tau, funk) ; 
end; 

Mtau: =vector(m+l) ; 

MtauCl] :=matrix( [CO ,0,0,0] , [0,0, 0,0] , [0,0, 0,0] , [0,0, 0,0]]) ; 
tau:= 0 ; 

M0 : =evalm(map(evaluera, integralen) ) ; 

for j from 1 to m 
do 

tau:=j*h/m; 

Mtau[j+1] :=evalm(map(evaluera,integralen)-M0) ; 
od; 

M : =Mtau [m+1] ; 

Minv: =evalm(inverse(M) ) ; 

ZZ : =evalm(Minv&*e_Ah.) ; 

WW : =evalm(transpose(e_Ali)&*ZZ+Minv) ; 
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WL: =submatrix(WW, [2,4] , [2,4] ) ; 

ZLU:=submatrix(ZZ, [2,4] , [2,4] ) ; 
ZLL:=transpose(submatrix(ZZ, [2,4] , [2,4] )) ; 

VR: =submatrix(WW, [2,4] , [1,3]); 

ZRU: =submatrix(ZZ, [2,4] , [1,3] ) ; 

ZRL: transpose (submatrix(ZZ, [1,3] , [2,4] )) ; 

xvel :=vector (n+1) ; 

xvel [1] : =evalm( (R [2] -R [l] ) /h) ; 

xvel [n+1] : =evalm( (R[n+1] -R[n] ) /b) ; 

Omega: =vector(n+l) ; 
for i from 2 to n 
do 

Omega [i] : =ZRL4*R[i - l] -WR&*R[i] +ZRU&*R[i+l] ; 
od; 

Omega [2] : =0mega[2] +ZLL£*xvel [1] ; 

Omega [n] : =Qmega[n] +ZLU&*xvel [n+1] ; 

DD : =vector(n+l) ; 

DD [2] : =WL ; 

for i from 3 to n 

do 

zd : =evalm(ZLL&* inverse (DD [i-1] ) ) ; 

DD [i] :=WL-zd&*ZLU; 

Omega [i] : =0mega[i] +zd&*Qmega[i-l] ; 
od; 

xvel [n] : =linsolve (DD [n] , Omega [n] ) ; 
for i from n-1 by -1 to 2 
do 

xvel [i] :=linsolve(DD [i] ,ZLU&*xvel[i+l]+Omega[i] ) ; 
od; 

x : =vector (n+1) ; 
for i from 0 to n 
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do 

x[i+l] :=matrix( [[R[i+1] [1,1]] , 

[xvel[i+l] [1,1]] , 

[R[i+1] [2,1]] , 

[xvel[i+l] [2,1]]]) ; 
od; 

plotlist : = [] ; 
for j from 1 to m 
do 

tau:=j*h/m; 

eAtau:=evalm(eAt(tau)) ; 

for i from 0 to n-1 
do 

entry : =evalm(eAtau&*(x [i+1] +Mtau [j+l]Jt*Minv 

&* (e_A [i+2] -x [i+1] ) ) ) ; 
plotlist : ={ [entry [1 , 1] , entry [3,1]], op(plotlist) } ; 
od; 
od; 

plot (plotlist, style=point) ; 
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